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ABSTRACT 

A  simple  derivation  of  Classman's  general  N  fast  Fourier  transform,  and 
corresponding  FORTRAN  program,  is  presented.  This  fast  Fourier  transform  is 
based  upon  a  representation  of  the  discrete  Fourier  transform  matrix  as  a 
product  of  sparse  matrices. 
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SIGNIFICANCE  AND  EXPLANATION 


The  discrete  Fourier  transform  is  the  basis  for  several  accurate 
techniques  for  the  numerical  solution  of  partial  differential  equations.  The 
fast  Fourier  transform,  an  alqorithm  which  allows  one  to  compute  rapidly  the 
discrete  FOurier  transform,  makes  these  techniques  computationally 
efficient.  This  paper  attempts  to  present  a  lucid  description  of  one  fast 
Fourier  transform,  the  fast  FOurier  transform  presented  by  Glassman. 

In  the  past  people  have  frequently  been  content  to  compute  rapidly  the 
discrete  FOurier  transform  of  vectors  whose  length  is  a  power  of  two. 
Glassman's  fast  Fourier  transform  allows  rapid  computation  of  the  discrete 
Fourier  transform  of  vectors  of  arbitrary  length. 


The  responsibility  for  the  wording  and  views  expressed  in  this  descriptive 
summary  lies  with  MFC,  and  not  with  the  author  of  this  report. 
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1 .  Introduction 


Let  the  N-vector  v 
i.e.,  the  components  v^ 


be  the  discrete  Fourier  transform  (DFT)  of  the  N-vector  u, 
of  v  are  computed  from  the  components  u^  of  u  hv  the  rule 


N 

=  I 

£=1 


(k-1  )U-1 ) 

U  (a) 

£  N 


for  k  =  1, 2,..., X 


where 


u)  =  exp{  -2tt/-T/N} 

N 

is  a  principle  N-th  root  of  unity.  It  is  easily  demonstrated  that  the  components  of  u 
can  be  recovered  from  the  components  of  v  by  the  rule 


1  N 
1  y 

n  v « 

k=1 


-(k-DU-1) 

7  U) 

k  N 


for  i  =  1,2,...,N 


The  N-point  DFT  matrix  WN  is  defined  to  be  the  matrix  of  order  N  whose  entry  in 
row  i,  column  j  is 

u(i-1)(j-1)  . 

Therefore  the  relations  between  u  and  v  presented  above  can  be  written  as 

v  =  W  u  and  u  =  -^  W„v 
N  N  N 

where  denotes  the  matrix  obtained  by  replacing  each  entry  of  WR  by  its  complex 
conjugate. 

A  fast  Fourier  transform  (FFT)  is  generally  considered  to  be  any  aloorithm  wh.ic'- 
rapidly  computes  the  DFT  of  a  given  vector.  One  of  the  most  popular  FFTs  was  presents-1  h- 
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Cooley  and  Tukey  [3]  in  1965.  Their  algorithm  commutes  the  DFT  of  an  N-vector  «sina 

N  .  <R  +  R2  +  ...  +  Rk) 

complex  operations,  where  one  operation  denotes  one  multiplicat ion  followed  by  one 
addition,  whenever  N  admits  the  representation 

N  =  P,R2  ...  Rk 

as  a  product  of  K  positive  integers  R1 , Rj# • • • » Since  the  publication  of  their 
article  numerous  authors  have  presented  other  FFTs,  each  requiring  approximately  the  sane 
number  of  complex  operations.  One  notable  exception  is  the  FFT  of  winoorad  [7] . 

In  this  paper  I  will  present  a  description  of  Glassman's  [5]  FFT.  This  description  of 
Glassman's  FFT  differs  from  one  presented  by  Drubin  f 4 }  only  in  the  definition  of  the 
tensor  product.  (However,  neither  Glassman  nor  Drubin  presented  a  FORTRAN  program  which 
computes  the  DFT  of  a  given  N-veotor . )  I  define  the  tensor  product  A  B  of  two 
matrices  A,  B  to  be  the  matrix  which,  when  partitioned  into  blocks  the  size  of  A,  has 
Abi  j  as  the  entry  in  block  row  i  and  block  column  j.  In  the  appendix  of  this  paoer  I 
have  presented  proofs  of  three  v*ell  known  properties  possessed  by  this  tensor  product. 

Glassman's  FFT  computes  the  DFT  of  an  N-vector  using  the  same  number  of  complex 
operations  as  the  Cboley-Tukey  FFT.  The  main  advantage  of  Glassman's  FFT  is  that  it  is 
easily  coded,  a  fact  which  should  be  compared  with  Singleton's  [6]  FFT.  The  main 
disadvantage  of  Glassman's  FFT  is  that  it  requires  an  N-vector  of  working  storage  to 
compute  the  DFT  of  an  N-vector.  I  will  show  how  one  can,  to  some  extent,  eliminate  this 
disadvantage. 

I  would  also  like  to  mention  that  de  Boor  [1J  has  recently  presented  an  FFT  that  is 
also  easily  described  and  coded. 


2.  Factorization  of  the  Discrete  Fourier  Transform  Matrix 


Consider  the  DFT  matrix  Wp^  where  P,Q  are  two  positive  integers. 
Partition  the  i-th  row  of  into  Q  groups  of  P  successive  entries, 

entries  in  the  q-th  group  are 

r  (i-1)(0+(q-1)P)  (i-1)(1+(q-1)P)  ( i-1 ) (P-1+(q-1 )P) • 

L  P  Q  '  Pp  PP 

Each  member  of  this  group  contains  the  common  factor 

( i-1 ) ( q-1 )P  (i-1)(a-1) 

<*)  =  tu  / 

PP  P 

therefore  the  q-th  group  admits  the  representation 

(i-U(q-l)  (P,P) 

“o 

where 

V(P.P)  -  r  ( i-1 ) ( 0  )  (i-1)(1)  ,  (i-1)(P-1)1 

yi  =  1“PQ  !  “PC  '  •  •  *f  “p Q  J 

denotes  the  first  p  entries  in  the  i-th  row  of  Wpp. 

Next,  partition  the  rows  of  into  P  groups  of  Q  successive  rows, 

the  p-th  group  are 


(  0  +(p-1)Q)(0)v(P,O)  .  <  0  +(p-1)Q)(C!-1)  (P,C>) 

B  #  T1+(P-1)Q . "B  ^  1+(p-1  )Q 

M(C-1  +  (P-1)C')(0)  (P,Q)  ,  (Q-1+(P-1)P)(Q-1)  (P,0) 

“b  V(P-1)Q . “o  V(P-i)p 

Observe  that  each  member  of  this  group  contains  the  term 

-  1 

Q  ' 

therefore  the  p-th  group  admits  the  representation 


Y(P.C) 

T 1+(p-1 )Q 


0 


{18 


V 


Here  the  matrix  in  square  brackets  is  a  block  diagonal  matrix,  each  block  a  1 
where  the  i-th  diagonal  block  is 


The 


The  rows  in 


x  P  matrix. 


<P.C> 

"  i+  (  n—  1  )r> 

and  Ip  is  the  identity  matrix  of  order  P. 

These  results  allow  us  to  prove  the  following 
Lemma :  The  DPT  matrix  Wp(5  admits  the  factorization 

Vc=f'p'q,{ip8  V 


(P,G>  r  ( i-1 )  ( 0 )  ( i- 1 )  { 1 )  (i-IHP-Di 

1 .  =  I  0) _  ;  <i) 

i  L  PC  PC  PC  J 

denotes  the  first  P  entries  in  the  i-th  row  of  wpp<  and 


1+(0)C 


1+(1)C 


1+(P-1)C- 


'C+(0)C 


Q+( 1 )Q 


'0+(P-l)C 


is  a  PC  x  0  block  matrix  with  1  x  p  blocks. 

Proof:  From  the  definition  of  we  find  that  the  p-th  group  of  C  successive  rows 


f<p,c>{ip  ®y 


y(P,0)  {Ip  8  W  }  , 

'V(p-i)cJ  w 


which  our  previous  computations  have  shown  to  be  the  p-th  group  of  0  successive  rows 
Wpg.  Since  p  was  arbitrary  we  therefore  conclude  that 

V  “  V  * 

The  matrix  F^P'^*  defined  in  the  above  lemma  has  several  interesting  limitino  ca 


F,P',»  -  W  and  F*1'^  =  I  . 

P  C 


in  particular 


These  observations  aid  us  in  the  proof  of  the  followina 
Theorem.  Let  N  admit  the  representation 


N  -  V2  •••  rk 


as  a  product  of  K  positive  integers  R., ,  R2, . . . ,  PK.  Then  WfI  admits  the  represe-.tat 

WN  -  F1F2  -  FK 

as  a  product  of  K  sparse  matrices  f.j  ,F2, . . .  ,fk  where 


Ft  =  I  „  0  F 

L  VV! 


( VRL+r**V 


(The  products  R1  ...  RL_1  for  L  =  1  and  RL+1  ...  By  for  L  =  K  are  defined  to  be  1 
Proof:  The  previous  lemma,  with  P  =  R1  and  ...  PK,  states  that 

WN  "  *\—*J  * 


Therefore  the  identity 


W  —  F  ...  F  ll  „  0  W„ 

N  1  L-1  VPL-1  V”Y 


holds  for  L  =  2.  Let  us  suppose  the  identity  holds  for  some  L  <  K.  The  previous  lemma 
with  p  =  Rl  and  Q  =  RL+1  ...  RR,  states  that 


W„  -  F 

V**rk 


(RL'RL+1”’RK> 


{ IR  8  WR  p  )  . 

rl  rl+i”'r* 


and  so 


IR  .  R  *  WR  R  “  PL{lR  R  8  WR  R  } 

1  "Vi  Y"RK  L  Rr,,RL  RL+1',,RK 


Consequently,  if  the  identity  holds  for  some  L  <  K  then  it  holds  for  L  +  1  too. 
Therefore  the  identity  must  hold  for  L  ■  K,  i.e. 


WN  ‘  F1F2  -•  FK 


where  we  have  noted  that 
I 


R.... 


a  w 


Vi  Rr  V”V 


(R„rD 

8  F  =  . 
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Let  C  >  1  be  a  divisor  of  B 
B  *  B/C 


u*  Ip  8F(C'A,u 


6. 

7. 

P. 

9.  endwhile  . 

With  the  exception  of  steps  6  and  8,  each  step  of  this  algorithm  can  be  directly 
implemented  in  FORTRAN.  Observe  that  step  6  admits  the  expansion 

6.1  C  *  2 

6.2  While  B  modulo  C  *  0  do 

6.3  C  =  C  +  1 

6.4  endwhile 

into  steps  that  can  be  directly  implemented  in  FORTRAN.  We  next  consider  the  expansion  of 
step  8. 

Let  the  product  RS  of  the  integers  R,S  be  a  divisor  of  N.  For  any  N-vector  w 
we  define  w*R*  to  be  the  FORTRAN  array  of  dimension  { R.N/R)  which  is  equivalent  to 
w,  and  w'"'°'  to  be  the  FORTRAN  array  of  dimension  (R,S,N/RS)  which  is  equivalent  to 


w.  This  definition  merely  implies  that 

<R) 


Wi.j  =  Wi+(j-1,R'  and 


Let 


(R.S) 

Wi,j,k  Wi+( j-1 )R+(k-1  )RS  ' 


v  -  iB  8F(C'A,u 


denote  the  result  of  the  computation  described  in  step  8.  As  shown  in  the  appendix  we  find 
that 

(B)  ( B)  _  ( C,  A  )T 


=  U  F 


or  equivalently  that 


,(p)  «  f  U(B)F(C,A) 

i.j  i-k  j,k 


for  i 


1,2 


B  and  j  =  1,2 


AC.  If  we  express  j  in  the  form 


j  =  ?A  +  (Dc  -  1>A  . 

with  1  <  j.  <  A  and  1  <  i  *  C,  then  the  nonzero  entries  in  tue  *•  -  r  *-  row  r* 
A  C 

are  the  numbers 

(  j-Df£-1) 

U) 

AC 

in  columns  k  *  t  +  (i  -  1  )r  f or  Jl  =  1,2,...,C.  We  thereforp  firv*  fchat 
A 

v(p>  _  ?  (p)  i,v,t(vi,Am-11 

Vi.jA+<  jc-1  )A  Ui,t  +  <iA-DC,JlAC 

or  equivalently  that 


v(E,A) 


C 


C 

I 

1  =  1 


(B.C)  <jA"1  +  <jC'1>RUi'1> 

u  to 

i,t,j,  AC 
A 


for  i  =  1,2, jA  =  1f2i***»A  and  =*  1,2,...,C.  Consecxuently,  step  °  admits 
expansion 


8.1 

For  jc 

=  1, 

8.2 

For 

^A 

=  1,2 . A 

8.3 

For 

x  — 

(B.A)  „  S  ( P,C)  (V 

i,jAjC  1-1  AC 

8.5 

Next 

i 

8.6  Next  jft 

8.7  Next  jc 

into  steps  that  can  be  directly  implemented  in  FORTRAN. 

Figure  1  presents  a  FORTRAN  version  of  Glassman's  FFT.  For  comparison  we 
de  Boor's  ft]  FFT  in  Figure  2.  I  have  found  that  Glassman's  FFT  runs  several  nercer.t 
faster  than  de  Boor's  FFT  on  the  University  of  Wisconsin's  UNI VAC  111°.  T^is  inrrea^® 
speed  is  probably  due  to  the  fact  that  the  loop  structure  used  in  rlass^ar's  rr?  nan  ~ 
efficiently  be  implemented  in  FORTRAN  than  the  loop  structu-e  used  in  4e  Four's  1 TT. 
increase  in  speed  would  therefore  vanish  if  one  were  to  hand  code  hot*'  rF7 '  s  us:'’  -i 
lanouage. 
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Appendix:  Some  Properties  of  the  Tensor  Product 


We  have  defined  the  tensor  product  A  8  P  of  two  matrices  A,B  as  the  matrix  which, 

when  partitioned  into  blocks  the  size  of  A,  has  Ab-  •  as  the  entry  in  block  row  i  and 

1 1  j 

block  column  j. 

Consider  now  any  N-vector  w.  If  R  is  a  divisor  of  N  we  define  w'  '  to  he  the 

FORTRAN  array  of  dimension  <R,N/R)  equivalent  to  w,  i.e. 

(R) 

”i,j  "  Wi+(j-1)R  - 

With  these  definitions  in  mind  let  us  now  prove  the  following 
Property  1:  Let  A,B  be  rectangular  matrices  where  A  is  a  R  x  c  matrix.  Then 

v  =  (A  8  B)u 


if  and  only  if 


(R) 


>  (C)=T 
Au  E 


Proof :  Let 


v  =  (A  8  B)u 


From  the  definition  of  the  tensor  product  A  8  B 

,<R)  _  v  ...  .  (C) 


we  observe,  for  each 
(C), 


i ,  that 


IK)  V  1C)  .ft  .  C) 

.  i  ■  2  “>i  *  »  b,  4u.  . 

j  i,j  '3  j  1,3  n 


The  sun  within  the  curly  brackets  is  easily  identified  as  the  i-th  column  of 


consequently  we  infer  that 


u(C>BT 


(R)  ,  (C)  T 

r  =  Au  B 


The  proof  of  the  converse  is  obtained  by  reversing  the  argument  presented  above.  • 

Carl  de  Boor  [2]  has  noted  that  this  property  allows  one  to  easily  compute 

v  «  (A  *  B)u 

given  u.  For  if  A  is  an  R  x  c  matrix  then 

(R)  »  (C>„T  J„,„  (C).T.T 

v  «  Au  B  «  |B(Au  )  J  f 

conseouently  programs  which  applv  A  and  B  to  vectors  can  easily  be  used  to  apply  A  8  P 
to  vectors.  This  property  also  allows  us  to  easily  prove  the  following 
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Property  2:  Let  the  products  A^Aj  an^  B1B2  defined.  Then 


ca1a2)  a  (p1b2)  =  (A1  a  a  b2)  . 


Proof:  Let  A^  be  an  x  matrix  for  k  =  1,2.  Observe  that  C1  =  Rj  because  the 

product  A-|*2  defined.  Let  u  be  an  arbitrary  vector  and  define 


w  =  i<A1A2>  8  (B1E2!  ' 


Property  1  implies  that 


(R,>  fC_ >  _  (Cl 

w  =  {A^Ju  2  <B1B2)T  =  A1(A2u  2  B2)B^  . 


If  we  define 


then  property  1  implies  that 


v  =  <a2  a  P2)u 


<V  «C2)  T 

v  =  A  u  B_  ,  and 

2  2 


<V  A  (V„T 

'  =  A1v  B1 


since  Ci  =  R2-  Using  property  1  once  more  we  find  that 

w  «  (A1  a  B1 )V  ,  and  so 

w  -  (a1  a  b1ha2  a  b2>u  . 

Consequently,  for  an  arbitrary  vector  u  *  have 


therefore 


{ ( A1  a2 )  a  (b1b2)}u  =  (A1  a  b1ha2  a  b2)u 


(a1a2)  a  (b1b2)  -  (A,  a  b^ia.,  a  b2>  . 


The  last  tensor  product  property  that  we  will  need  is  described  as  follows. 
Property  3:  For  arbitrary  matrices  A1(  A2  and  Aj 


(a2  a  a3>  »  <a1  a  a2 )  a  a3  . 


Proof:  I*«t  A^  be  an  x  ck  matrix  for  k  «  1,2,  and  3.  Let  e^B'  be  the  B-vector 

obtained  by  replacing  the  i-th  component  of  the  zero  B-vector  by  1.  Let  a(k).  denote  the 

i ,  j 

entry  of  Ak  in  row  i  and  column  j .  Observe  that 


(F^l  <C,  )T 


a  =  l  a'k’e.  e.  k  for  V  =  1,2.3 

\  x.J  i  3 


Consequently 


IP  )  ( C  IT  IP,)  "V" 


a,  a  <a2  a  a3>  .  i  .ixw«[.r,,.r,,T:  •  ((•r2*?”’  •  •% ?e, 3 


and 


(A 


(  a  #j)  a  A3  ■  S  *i]jak!!lani!n^ei  *j  '  ^  8  ^ek  ^  ei 


(R.)  (C,)T  (R,>  (C  )T  .  ( R, )  (C  )T 

I  I  T  f  *  ^  *  7  1  -  J  J 


e  e 
k  m  n 


From  the  easily  verified  identity 


(R.)  <C  >T  (Rj)  (C,)T,  r  (P,)  <C3)Ti 

•j  1  •  (K  ei 


]»^3en3  H- 


(RJ  ( C  )T  #  <R,>  (C  )T  r  <M  <CJT 


([«i  ej  1  *  ^6k  ®i 


1 )  ®  [e  e 

■* '  L  m  n 


A1  *  (A2  *  A3>  =  (A1  *  A2*  *  A3 


we  deduce  that 
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